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Abstract. We analytically investigate the nonlinear response of a damped doubly 
clamped nanomechanical beam under static longitudinal compression which is excited 
to transverse vibrations. Starting from a continuous elasticity model for the beam, we 
consider the dynamics of the beam close to the Euler buckling instability. There, the 
fundamental transverse mode dominates and a quantum mechanical time-dependent 
effective single particle Hamiltonian for its amplitude can be derived. In addition, 
we include the influence of a dissipative Ohmic or super-Ohmic environment. In the 
rotating frame, a Markovian master equation is derived which includes also the effect of 
the time-dependent driving in a non-trivial way. The quasienergies of the pure system 
show multiple avoided level crossings corresponding to multiphonon transitions in the 
resonator. Around the resonances, the master equation is solved analytically using 
Van Vleck perturbation theory. Their lineshapes are calculated resulting in simple 
expressions. We find the general solution for the multiple multiphonon resonances 
and, most interestingly, a bath-induced transition from a resonant to an antiresonant 
behavior of the nonlinear response. 
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1. Introduction 

The experimental realization of nanoscale resonators which show quantum mechanical 
behavior [TJ 12 El IH E] is currently on the schedule of several research groups worldwide 
and poses a rather non-trivial task. Important key experiments on the way to this 
goal have already been reported in the literature pi7|IH|I^ITTl | ITT | IT2 l lT ^ rfI | ITB l lTn 1 
IT71 HH1 EH 120] an d are also reviewed in this Focus Issue. Most techniques to reveal 
the quantum behavior so far address the linear response in form of the amplitude of 
the transverse vibrations of the nanobeam around its eigenfrequency. The goal is to 
excite only a few energy quanta in a resonator held at low temperature. To measure 
the response, the ultimate goal of the experiments is to increase the resolution of the 
position measurement to the quantum limit jTTJ [T7J UHl I2H 122] • As the response of a 
damped linear quantum oscillator has the same simple Lorentzian shape as the one of a 
damped linear classical oscillator 23 , a unique identification of the "quantumness" of 
a nanoresonator in the linear regime can sometimes be difficult. 

One possible alternative is to study the nonlinear response of the nanoresonator 
which has been excited to its nonlinear regime. A macroscopic beam which is clamped 
at its ends and which is strongly excited to transverse vibrations displays the properties 
of the Duffing oscillator being a simple damped driven oscillator with a (cubic) nonlinear 
restoring force j21j. Its nonlinear response displays rich physical properties including 
a driving induced bistability, hysteresis, harmonic mixing and chaos |2H 1221 12H] • The 
nonlinear response of (still classical) nanoscale resonators has been measured in recent 
experiments [HI GUI EJ EH]- In the range of weak excitations, the standard linear response 
arises while for increasing driving, the characteristic response curve of a classical Duffing 
oscillator has been identified. 

No signatures of a quantum behavior in the nonlinear response of realized 
nanobeams have been reported up to present. One reason is that a nanomechanical 
resonator is exposed to a variety of intrinsic as well as extrinsic damping mechanisms 
depending on the details of the fabrication procedure, the experimental conditions and 
the used materials |2lH |23 123 121] ■ Possible extrinsic mechanisms include clamping 
losses due to the strain at the connections to the support structure, heating, coupling 
to higher vibrational modes, friction due to the surrounding gas, nonlinear effects, 
thermoelastic losses due to propagating acoustic waves, surface roughness, extrinsic noise 
sources, dislocations, and other material-dependent properties. An important internal 
mechanism is the interaction with localized crystal defects. Controlling this variety of 
damping sources is one of the major tasks to be solved to reveal quantum mechanical 
features. Recent measurement show that in the so far realized devices based on silicon 
and diamond structures, damping has been rather strong at low frequencies [23 I2H] 
indicating even sub-Ohmic-type damping j2H| which would make it difficult to observe 
quantum effects at all. However, using freely suspended carbon nanotubes [3U1 EI] 
instead could reduce damping at low frequencies due to the more regular structure of 
the long molecules which can be produced in a very clean manner. Further experimental 
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work is required to clarify this point and to optimize the experimental conditions. 

Nanoscale nonlinear resonators in the quantum regime have been investigated 
theoretically starting from microscopic models based on elasticity theory for the beam 
[321 EH1 E3! • Carr, Lawrence and Wybourne have considered an elastic bar under static 
longitudinal compression beyond the Euler instability leading to two stable equilibrium 
positions around which the transverse vibrations of the beam occur. It turned out that 
quantum tunneling between the two minima is in principle possible in silicon beams 
and carbon nanotubes. However, the strain has to be controlled with extreme accuracy 
and the quantum fluctuations in position are of the order of 0.1 A. The detection of 
such small lengths certainly is challenging. However, a possible method to increase 
the resolution could be the use of the phenomenon of stochastic resonance [33] for a 
coherent signal amplification of the nonlinear response of nanomechanical resonators in 
their bistable regime [3E] • Werner and Zwerger [SSj have studied a similar setup close to 
the Euler buckling instability which occurs at a critical strain e c . There, the frequency 
of the fundamental mode vanishes and quartic terms in the Lagrangian have to be 
taken into account. An effective Hamiltonian has been derived for the amplitude of 
the fundamental mode being the dynamical variable which moves in an anharmonic 
potential. Depending on the strain e being below (e < e c ) or above (e > e c ) the 
critical value, a monostable or bistable situation can be created. The conditions for 
macroscopic quantum tunneling to occur have been estimated for the bistable case. In 
order to measure single-phonon transitions in a nanoresonator, it has been proposed to 
use its anharmonicity together with a second nanoresonator acting as a transducer for 
the phonon number in the first one [37]. In this way, the measured signal being the 
induced current is directly proportional to the position of the read-out oscillator. 

In Ref. [SI], we have considered a similar setup but restricted to the statically 
monostable case below the Euler instability, i.e., for e < e c . In addition, we have allowed 
for a time-dependent periodic driving force F(t) such that an effective monostable 
quantum Duffing oscillator arises. Possible origin of the driving can be the magneto- 
motive force when an ac current is applied and the beam is placed in a transverse 
magnetic field. Moreover, a (weak) influence of the environment has been modeled 
phenomenologically by a simple Ohmic harmonic bath. The nonlinear response has 
been determined numerically from solving a Born-Markovian master equation for the 
reduced density operator of the system after the bath has been traced out. We have 
identified discrete multiphonon transitions as well as macroscopic quantum tunneling 
of the fundamental mode amplitude between the two stable states in the driving 
induced bistability. Moreover, a peculiar multiphonon antiresonant behavior has been 
found in the numerical results for the damped system [3H]. The discrete multiphonon 
(anti-)resonances are a typical signature of quantum mechanical behavior |3i] l3*5] and are 
absent in the corresponding classical model of the standard Duffing oscillator [231EE1EII], 
also when thermal fluctuations are included [35] . 

While we have approached the problem in Refs. |3i] EH] by numerical means, we 
present in this work a complete analytical investigation of the dynamics of the quantum 
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Duffing oscillator. We intend to elucidate the mechanism behind the reported [SB] 
bath-induced transition from the resonant to the antiresonant nonlinear response of 
the nanobeam. This is achieved by solving a Born-Markovian master equation for the 
reduced density operator in the rotating frame. Within the rotating wave approximation 
(RWA), a simplified system Hamiltonian follows whose eigenstates are the quasienergy 
states. The corresponding quasienergies show avoided level crossings when the driving 
frequency is varied. They correspond to multiphonon transitions occurring in the 
resonator. Moreover, we include the dissipative influence of an environment and find 
that the dynamics around the avoided quasienergy level crossings is well described 
by a simplified master equation involving only a few quasienergy states. Around the 
anticrossings, we find resonant as well as antiresonant nonlinear responses depending on 
the damping strength. The underlying mechanism is worked out in the perturbative 
regime of weak nonlinearity, weak driving and weak damping. There, Van Vleck 
perturbation theory allows to obtain the quasienergies and the quasienergy states 
analytically. The master equation can then be solved in the stationary limit and 
subsequently, the line shapes of the resonant as well as the antiresonant nonlinear 
response can be calculated. 

The problem of a driven quantum oscillator with a quartic nonlinearity has been 
investigated theoretically in earlier works in various contexts. In the context of 
the radiative excitation of polyatomic molecules [UI|, Larsen and Bloembergen have 
calculated the wave-functions for the coherent multiphoton Rabi precession between 
two discrete levels for a collisionless model. More recently, also Dykman and Fistul [H] 
have considered the bare nonlinear Hamiltonian under the rotating wave approximation. 
Drummond and Walls [12] have investigated a similar system occurring for the case of 
a coherently driven dispersive cavity including a cubic nonlinearity. Photon bunching 
and antibunching have been predicted upon solving the corresponding Fokker-Planck 
equation. Vogel and Risken [U"| have calculated the tunneling rates for the Drummond- 
Walls model by use of continued fraction methods. Dmitriev, D'yakonov and Ioffe 
[4*4] have calculated the tunneling and thermal transition rates for the case when 
the associated times are large. Dykman and Smelyanskii [f"""J have calculated the 
probability of transitions between the stable states in a quasi-classical approximation in 
the thermally activated regime. Recently, the role of the detector (in this case, a photon 
detector) has been studied for the quantum Duffing oscillator in the chaotic regime 
[IH]. The power spectra of the detected photons carry information on the underlying 
dynamics of the nonlinear oscillator and can be used to distinguish its different modes. 
However, the line-shape of the multi-phonon resonance which is the central object for 
studying the nonlinear response remained unaddressed so far. In addition, we start 
from a microscopic Hamiltonian for the bath and present a fully analytical treatment 
of system and environment in the deep quantum regime of weak coupling. 

Our paper is structured such that we introduce the elasticity model for the doubly 
clamped nanobeam, derive the effective quartic Hamiltonian and discuss the model for 
damping in Section ""0 Then, we discuss the coherent dynamics of the pure system in 
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terms of the RWA and the Van Vleck perturbation method in Section |H1 The dissipative 
dynamics is studied in Section while the observables are defined in Section |SJ The 
solution for the line shapes are given in Section El before the final conclusions are drawn 
in Section [7| 



2. Model for the driven nanoresonator 



We consider a freely suspended nanomechanical beam of total length L and mass 
density a = m/L which is clamped at both ends {doubly clamped boundary conditions) 
and which is characterized by its bending rigidity \x = YI being the product of 
Young's elasticity modulus Y and the moment of inertia /. In addition, we allow 
for a mechanical force F > which compresses the beam in longitudinal direction. 
Moreover, the beam is excited to transverse vibrations by a time-dependent driving 
field F{t) = fcos(u cx t). In a classical description, the transverse deflection (fi(s,t) 
characterizes the beam completely, where < s < L. Then, the Lagrangian of the 
vibrating beam follows from elasticity theory as (SB] 

L r a. 2 /i 0" 2 



C(4>, 4>,t) = / ds 







V 7 ! - 0' 2 - l) + F(t)<f> 



T 2 1-0' 2 

Before we study the dynamics of the driven beam, we consider first the undriven system 
with Fit) = 0. 

For the case of small deflections |0'(s)| <C 1, the Lagrangian can be linearized 
and the time-dependent Euler-Lagrange equations can be solved by the eigenfunctions 
(f)(s,t) = ^2 n (fin(s,t) = J2 n A n (t)g n (s), where g„{s) are the normal modes which follow 
as the solution of the characteristic equation. For the doubly clamped nano-beam, we 
have 0(0) = 4>{L) = and 0'(O) = 4>'(L) = 0. However, it turns out that this situation is 
closely related to the simpler case that the nano-beam is also fixed at both ends but its 
ends can move such that the bending moments at the ends vanish, i.e., 0(0) = <p{L) = 
and 0"(O) = 4>"{L) = [free boundary conditions). For the case of free boundary 
conditions, the characteristic equation yields the normal modes g^ ee (s) = sm(nirs/L) 
and the corresponding frequency of the n— th mode follows as 

free _ ( ^71 / L) 2 - F \ 1/2 UTX 

Un -{ ; ) T* (2) 

At the critical force F c = ^{tt /L) 2 , the fundamental frequency a;f ee (Fo — > F c ) vanishes 
as \fe, where e = (F c — F )/F c is the distance to the critical force, and the well-known 
Euler instability occurs. 

For the case of doubly clamped boundary conditions, the characteristic equation 
yields a transcendental equation for the normal modes which cannot be solved 
analytically. However, close to the Euler instability Fq —> F c , the situation simplifies 
again. After expanding, one finds for the fundamental frequency Ui(F — > F c ) = i/ecuo, 
with the frequency scale ujq = (4/v^) \/pT/a(Tc/L) 2 . Approaching the Euler instability, 
the frequencies of the higher modes u n >2 remain finite, while the fundamental frequency 
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uj\{Fq — > F c ) vanishes again like \fe. Hence, the dynamics at low energies close to the 
Euler instability will be dominated by the fundamental mode alone which simplifies the 
treatment of the nonlinear case, see below. The fundamental mode gi{s) can also be 
expanded close to the Euler instability and one obtains in zero-th order in e 

^(s) ~ sin 2 (^) . (3) 
2.1. Effective single-particle Hamiltonian 

Since the fundamental mode vanishes when Fq — > F c , one has to include the 
contributions beyond the quadratic terms oc 0' 2 , <fi" 2 of the transverse deflections in the 
Lagrangian. The next higher order is quartic and yields terms oc /4 ,0 /2 0" 2 . Inserting 
again the normal mode expansion in the Lagrangian generates self-coupled modes J2k 
as well as couplings terms J2k i A\A 2 between the modes. This interacting field-theoretic 
problem cannot be solved any longer. However, since the normal mode dominates the 
dynamics at low energies closed to the Euler instability, one can neglect the higher modes 
in this regime. Hence, we choose the ansatz <p(s,t) = A\(t)gi(s) in the regime Fq — > F c 
and restrict the discussion in the rest of our work to this regime. The so-far classical 
field theory can be quantized by introducing the canonically conjugate momentum 
V = —ihd/dAi and the time-dependent driving force can straightforwardly be included. 
Note that when the driving frequency is close to the fundamental frequency of the beam, 
the fundamental mode will dominate also in absence of a static longitudinal compression 
force. However, a compression force helps to enhance the nonlinear effects which are 
in the focus of this work. After all, an effective quantum mechanical time-dependent 
Hamiltonian results which describes the dynamics of a single quantum particle with 
"coordinate" X = Ai in a time-dependent anharmonic potential. It reads 

with the effective mass m e g = 3aL/8 and the nonlinearity parameter a = (tt/L) 4 F c L(1 + 
3e). The classical analogous system is the Duffing oscillator [21] (when, in addition, 
damping is included, see below). It shows a rich variety of features including regular 
and chaotic motion. In this work, we focus on the parameter regime where only regular 
motion occurs. For weak driving strengths, the response as a function of the driving 
frequency u; ex has the well-known form of the harmonic oscillator with the maximum at 
u ex = u)\. For increasing driving strength, the resonance grows and bends away from 
the o; ex = co>i-axis towards larger frequencies (since a > 0). The locus of the maximal 
amplitudes is often called the backbone curve [24J. The corresponding nonlinear response 
of the quantum system shows clear signatures of sharp multi-phonon resonances whose 
line shapes will be calculated below. 
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2.2. Phenomenological model for damping 

In our approach, we do not intend to focus on the role of the microscopic damping 
mechanisms as this depends on the details of the experimental device. Instead, 
we introduce damping phenomenologically in the standard way [22] by coupling the 
resonator Hamiltonian Eq. (j3J) to a bath of harmonic oscillators described by the 
standard Hamiltonian 

"- „:v)'- < 5 > 

with the spectral density 

c 2 

■7H = ? Yl - = m eS rr.ul-'u'e- u ' u ° , (6) 

with damping constant 7 S and cut-off frequency uj c . Our results discussed below are 
valid for an Ohmic (s = = 7) as well as for super-Ohmic (s > 1) baths. Sub- 
Ohmic baths will not be considered here since the weak-coupling assumption which 
allows the Markov approximation does not hold any longer. Formally, the coefficients 
in the master equation would diverge in the sub-Ohmic case, see Eq. (J2H|) below. The 
total Hamiltonian is H tot (t) = H(t) + Hb- 

To proceed, we scale H tot {t) to dimensionless quantities such that the energies are 
in units of Hui while the lengths are scaled in units of xq = \/ ~^~ ■ Put differently, 
we formally set m c s — K — iO\ — 1. The nonlinearity parameter a is scaled in units of 
a = tkjji/ , while the driving amplitudes are given in units of f = hLo 1 /x . Moreover, 
we scale temperature in units of T = hui/k B while the damping strengths are measured 
with respect to u 1 . 



3. Coherent dynamics and rotating wave approximation (RWA) 

Let us first consider the resonator dynamics without coupling to the bath. For 
convenience, we switch to a representation in terms of creation and annihilation 
operators a and a), such that X = xo(a + a))/\/2. Moreover, it is convenient to 
switch to the rotating frame by formally performing the canonical transformation 
R = exp [— iu; ex aW] . We are interested in the nonlinear response of the resonator around 
its fundamental frequency, i.e., for u ex w ui, and will not consider the response at higher 
harmonics. We further assume that the driving amplitude / is not too large such that 
the nonlinear effects are small enough in order not to enter the chaotic regime. This 
suggests to use a rotating wave approximation (RWA) of the full system Hamiltonian 
H(t) in Eq. (@J) as the fast oscillating terms will be negligible around the fundamental 
frequency for weak enough driving. By eliminating all the fast oscillating terms from the 
transformed Hamiltonian, one obtains the Schrodinger equation in the rotating frame 
H\<fi a ) = e a \(f> a ) with the Hamiltonian in the RWA 

H = ojfi+ V -h{n + 1) + / (a + a}) . (7) 
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Here, we have introduced the detuning Co = uj\ — cj cx , the nonlinearity parameter 
v = 3ha I (4ujf) , / = f^tvjJi)^ 1 ! 2 and n = a)a. In the static frame, an orthonormal 
basis (at equal times) follows as 

\<p a (t)) = e-^ at \4> a ) . (8) 

The Hamiltonian (J7j) has been studied in Refs. [1UJ E]. The quasi-energy levels for a 
given number N of phonons are pairwise degenerate, £jy-n = e n for n < N, vanishing 
/ — > and u) = —u(N + l)/2. For a finite driving strength / > 0, the exact crossings 
turn then into avoided crossings which is a signature of multiphonon transitions [SUE]. 
A typical quasienergy spectrum is shown in Fig. ^ for the parameters v = 10 -3 and 
/ = 10~ 4 . The dashed vertical lines indicate the multiple avoided level crossings which 
occur all for the same driving frequency. For |e| = \2f /[u(N + 1)]| <C 1, each pair of 
degenerate levels interacts only weakly with the other levels, and act effectively like a 
two- level Rabi system [12]. The Rabi frequency is related to the minimal splitting of 
the levels and is calculated perturbatively with e as a small parameter in the following 
Section. 



3.1. Van Vleck perturbation theory 



Let us therefore consider the multiphonon resonance at uj = —v(N + l)/2. In addition, 
we are interested in the response around the resonance and therefore introduce the small 
deviation A. We formally rewrite H as 



H 



u(N + l) 



-a 



A)n + — — n + e(a + a 1 ) 



N + l 

Let us then first discuss the dynamics at resonance (A = 0) 
unperturbed part H and the perturbation eV according to 



(9) 

We divide it in the 



n+l 
-n + — — T n 



V 



[a + a'] 



(10) 



+ 1 

respectively. The unperturbed Hamiltonian is diagonal and near the resonance 
its spectrum is divided in well separated groups of nearly degenerate quasienergy 
eigenvalues. An appropriate perturbative method to diagonalize this type of 
Hamiltonian is the Van Vleck perturbation theory |S01 EH]- It defines a unitary 
transformation yielding the Hamiltonian H in an effective block diagonal form. The 
effective Hamiltonian has the same eigenvalues as the original one, with the quasi- 
degenerate eigenvalues in a common block. The effective Hamiltonian can be written 
as 



H' 



AS 



He 



-iS 



In our case, each block is a two by two matrix corresponding to a subspace formed by 
a couple of quasienergy states forming an anticrossing. Let us consider the effective 
Hamiltonian H' n corresponding to the involved levels |n) and \N — n), being eigenstates 
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of the harmonic oscillator. The degeneracy in the corresponding block is lifted at order 
iV — In in Van Vleck perturbation theory. The block Hamiltonian then reads 

fr,_( |n(n-JV) e N - 2n C 1%N . 2n \ 

n \ e N - 2n C l2 , N . 2n fn(n - N) )> [iZ) 

where 

c u N - 2n = (n + i)^ . (13) 

This is the lowest order of the perturbed Hamiltonian which allows to calculate the 
corresponding zero-th order eigenstates. By diagonalizing H' n in Eq. (|12j). one finds the 
minimal splitting for the N— phonon transition as 

&N,n = \2 £N 2U Ci2,N-2n\ 



= 2f (v) K -^ 1 ^EE«L. (14) 

V v ) y/nl(N-2n- l)! 2 

For the case away from resonance, we consider a detuning A = e N 5. Within the Van 
Vleck technique, only the zero-th block is influenced according to 



„/ s n C 12 ,n \ n 



the other blocks given in Eq. (JT2j) for n ^ are not influenced by this higher-order 
correction. The eigenvectors for the Hamiltonian H at zero-th order are obtained 
by diagonalizing H' n in Eq. (|12|). One finds \<f> n ) = \n) for n > N + 1 or |0 n ) = 
(\n) + \N-n))/V2&nd \<p N - n ) = {\n)-\N-n))/V2ior < n < N/2 and \<j> N/2 ) = \N/2) 
if N is even. Moreover, 

|0 O ) = cos^|0)-sin^|iV), 

|0^)= sin^|0) + cos^|V), (16) 
where we have introduced the angle 9 via tan# = —2Qn^/[i , (N + 1)NA]. 



4. Dissipative dynamics in presence of the bath 

Having discussed the coherent dynamics, we include now the influence of the harmonic 
bath coupled to the driven system. We therefore assume that the coupling is weak 
enough such that the standard Markovian master equation 

j t Q = -i[H(t) 1 Q]+C Q (17) 

for the reduced density operator p(t) can be applied. The influence of the bath enters 
in the superoperator 

Cq = -[X, [P(t), £?]+] - [X, [Q(t), g)) (18) 
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(19) 



(20) 



7 (r) = 



J{u) 



irm eS j u 
1 f°° 

K(t) = - / J{uj) coth 

JO 



COS LOT , 

2k B T 



COS U1T , 



(21) 



where T is the environment temperature. Moreover, U(t,t') = Texp(i f*H(t)dt) is the 
propagator with the time order operator T. Next, we project the density matrix on the 
orthonormal set \(p a (t)) = exp [— iuj eyL a)at] \(f> a ), such that the matrix elements read 

Q«p(t) = (0a|e^°S(*)e- i ^ otat |0 / 9> • 
Performing the derivative one obtains 



Q a p(t) = -i(<f>, 



He 



ra^at 



i + [H(t), g] + iCg + gi^ 



-iu) cy -(V at I 



~ -i{e a -ep) eaP {t) + (Mt)\e tuJ ^ at Cge-^ at 
For the dissipative term, we need to compute 



\Mt))- (22) 



X aP (t) = (0 a \e luJc * a ' at Xe- luJ ™ atat \(f)p) 

= (e-^X a ^ +1 + e^X^) , (23) 

with X a p >+ i = Xp i_ = x o (0 a |a|0 / 3)/v / 2 being the matrix element of the destruction 
operator in the rotating frame. We need, moreover, 

POO 

Qaf3(t) = / drK{r){(t> a \e^ at U\t - r,t)XU(t - r^e'^^p) 
Jo 

poo 

= / rfrir(r)e- i(£ "- £ ' 3)r (0 a |e^ xata( *- r) A'e-^ cxata( *- r) |0 /3 ) (24) 
Jo 

' poo 

/ dTK(r)e- i{£a ~ £ p-^ )T 
Jo 



= e 



-iu> cx t 



a/3,+1 



+ e 



drK(r)e~ 



i(e a — £ / g+w e x)r 



X, 



a/3-1 



(25) 



In an analogous way, we have 

1 f°° + 
P a(j (t) = - dT-f(r)e-^ a - £ ^ T ((j) a \e" a{t ' T) Ve- 

2 Jo 



10/3) 



poo 

/ dr7(r)< 
Jo 



-i(ea-£/3)r 
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x [ (0 Q |e~ ta(i - T) [±z j +H{t),X]e- l ^ a]a{t ~ T) \(t)p) 



d . 

-j t { " 



| e iw ox at a(t-r) e ~iuj c ^ a(t-r) | 

Sa-ep- ij^j J™ dr^e-^-^X^t - r) 



(26) 



Here, we have defined the time derivative ± 4 where the positive (negative) sign belongs 
to the left (right) direction. In the second line, we have used the canonical relation 
'P/ m eS — —i[X,H(t)]. We can now compute the matrix elements of the operators 
involved in the dissipative part in Eq. (|22j) and find for the terms in Eq. (|18|) 

( p + Q)ap = (e'^Nafr-iXa^+i + e^N^+xX^x) (27) 

and 

(P ~ Q) a p = ~ {e-^Nfa+xX^+x + e^N^.xX^x) . (28) 
Here, N a p ; ±x are defined as 

N Ql 3,±x = N(e a -e p ± u ex ) N(e) = J(\e\)[n th (\s\) + 9{-e)\ , (29) 

in terms of the bath density of states the bosonic thermal occupation number 

1 



n th(e) , } 



COthl 2feT 1 - 1 



(30) 



and the Heaviside function 9(x). 

Eq. ()29|) illustrates why we have to restrict to (super-) Ohmic baths, since N(e) 
would diverge for s < 1 at low energies. During the calculation, the r-integration in the 
double integrals in Eqs. (|2*Sjl and (|2~5j) has been evaluated by using the representation 
J °° drexp^iur) = tx5{uj) + iV p (l/w), where V p denotes the principal part. The 
contributions of the principal part result in quasienergy shifts of the order of 7 S which 
are the so-called Lamb shifts. As usual, these have also been neglected here. 

The ingredients can now be put together to obtain the Markovian master equation 
in the static frame as 



E E 



e 



i(n+n )u>ext \^]Sf aa ,_ n -f jV 'pp', n ')X aa i >n Pa/3' X pi p^ n i 



a'P> n,n'=±l 

— Npi a i- n iX a p^ n Xpi a iyp a ip — IS ' a i pi <n > p a pi X pi a i )fl i X a i p jTl ] . (31) 

Next, we perform a 'moderate rotating wave approximation' consisting in averaging the 
time-dependent terms in the bath part over the driving period T Wex = 2n/u) ex . This is 
consistent with the assumption of weak coupling which assumes that dissipative effects 
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on the dynamics are noticeable only on a time scale much larger than T UeK . Under this 
approximation, the master equation becomes 

Qap{t) — Sa/3,a'0' Qa'P 1 jt) 
a'/3' 

= [— i(e a — Sp)8 aa ib~ppi + CaP,a'P')]Qa'P'{t) (32) 

a'P' 

with the dissipative transition rates 

£ a p,a>P' = {Naa'-ri + Nppi - n )X aa i )n Xpip- n 

n=±l 

— 8aa' ^a"P'-nXp' a ",-nX a "p,n 
a";n=±l 

~ fipP' T P" 'a' ,-nX a p" - n Xpii a > >n . (33) 
p" ;n =±l 

It is instructive to compare this master equation to the one in Ref. given in terms 
of the full Floquet quasienergy states. The key difference here is that the density matrix 
is projected onto the approximate eigenvectors exp (— iu ex a^at)\(p a ) rather than onto the 
exact Floquet solutions. As a consequence of the RWA, the sums in Eq. only include 
the n = ±1 terms indicating that only one-step transitions are possible where n = — 1 
refers to emission and n — +1 to absorption. Being consistent with the RWA, we can 
assume that \u\, \f\, \u cx — uji\ <^ uj\ which yields to \e a — ep\ u ex . Hence, N a p }+ i is the 
product of the bath density of states and the bosonic occupation number at temperature 
T. This corresponds to the thermally activated absorption of a phonon from the bath. 
On the other hand, N a p-i given in Eq. (|2!?|) contains the temperature-independent term 
.. + J(u; cx ) describing spontaneous emission. 

5. Observable for the nonlinear response 

Assuming an ergodic dynamics of the full system, or equivalently that there is just 
one eigenvector of the superoperator S in Eq. (}32|) . corresponding to a vanishing 
eigenvalue, and that all the other eigenvalues have negative real part, the asymptotic 
solution of Eq. (J3~2~|) is 

lim g(t) = q 00 . (34) 

t^oo 

To simplify notation, we omit in the following the superscript oo but only refer to the 
stationary state g a p = g^. 

We are interested in the mean value of the position operator in the stationary state 
according to 

(X) t = tr(gX) = 6apXp a (t) . (35) 

aP 

Using Eq. (J2l?|) yields (X) = A cos (uj ex t + <p), with the oscillation amplitude 

A = 2 1 ^ e»P X p a:+1 1 , (36) 

aP 
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Re 







and the phase shift 

ip = tt6 J —Re Q a pXp a ,+i J + arctan |= ^ > (37) 

V La/3 J / 

with # being the Heaviside function. 

6. Analytical solution for the lineshape of the multiphonon resonance in 
the perturbative regime 

When the driving frequency u> ex is varied, the amplitude A shows characteristic multi- 
phonon resonances at those values for which the quasienergy levels form avoided level 
crossings [33]. While in Ref. [33] these resonances have been studied numerically, it is 
the central result of this work to calculate their line shape analytically by solving the 
corresponding master equation in the Van Vleck perturbative regime. 

Within the limit of validity of the RWA, i.e., \u\, \f\, \u cx — u>i\ <C u>i, we have 
k« — £p\ "C uj ex . In the regime of low temperature k B T <C u cx , it follows from Eq. 
(|29|) that Nafi-i ~ J(uJ e x) and N a p t x ~ entering in the transition rates in Eq. (|3*3*|) . 
This approximation corresponds to consider spontaneous emission only and yields the 
dissipative transition rates 



13, a' 13' 



2 



13/3' 



a"P'A a iip 



~&PP' Ap"aAf3" a ' J 

P" J 



(3f 



Here, we have defined A a p = (<fi a \a\<f)p). Note that it is consistent with the previous 
approximation to set uj cx /uji « 1. Hence, all the following results are valid for Ohmic 
as well as super-Ohmic baths. 

In the following we will use this simplified transition rates to solve the master 
equation near the multiple multiphonon resonances. The transition between the 
groundstate and the iV-phonon state is the narrowest. Hence, it will be affected first 
when a finite coupling to the bath is considered. In particular, it is interesting to 
consider the case when the damping constant 7 S is larger than the minimal splitting 
Qnq between the two quasienergy states but smaller than all the minimal splittings of 
the other, i.e., f2/v < 7s "C for n > 1. In this case, we can assume a partial 
secular approximation: We set all the off-diagonal elements to zero except for q n and 
Qno — Qon- m this regime the stationary solutions are determined by the conditions 

= £-aa,PpQpp + C-aaflN 2 Re(^ov) , 

p 

(39) 

a 

For very weak damping, i.e., when 7 S is smaller than all minimal splittings 
(7 S <C fijvn), the off-diagonal elements of the density matrix are negligibly small and 
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can be set to zero. Within this approximation, the stationary solution for the density 
matrix is determined by the simple kinetic equation 

= y.Cgg^pQpp. (40) 

P 

In this regime, a very simple physical picture arises. The bath causes transitions between 
different quasienergy states, but here, the transition rates are independent from the 
quasienergies. It is instructive to express the quasienergy solutions in terms of the 
harmonic oscillator (HO) solutions as \4> a ) = J2 n c an\ n ) with some coefficients c an . The 
transition rates between two quasienergy states then read 

C aa ,pp = 'J s \((j) a \a\(f) f3 )\ 2 = 7 s ^(n + l)\c an \ 2 \cf3n+i\ 2 ■ (41) 

II 

This formula illustrates simple selection rules in this low-temperature regime: Only 
those components of the two different quasienergy states contribute to the transition 
rate whose energy differs by one energy quantum {n «-> n + 1). 

6.1. One-phonon resonance vs. antiresonance 

Before we consider the general multiphonon case, we first elaborate on the one-phonon 
resonance. This, in particular, allows to make the connection to the standard linear 
response of a driven damped harmonic oscillator which is resonant at the frequency 
Ui + v. We will illustrate the mechanism how this resonant behavior is turned into an 
antiresonant behavior when the damping is reduced (and the driving amplitude / is 
kept fixed). 

The corresponding effective Hamiltonian H' Q follows from Eq. (|15jl and is readily 
diagonalized by the quasienergy states |0o) and |0i) which are of zero-th order in e 
and which are given in Eq. ()16j) . The master equation ()39|) can be straightforwardly 
solved in terms of the rates C a p^ a 'i3' for which one needs the ingredients A o = —A\\ = 
sin(6>/2)cos(6>/2),v4 01 = cos 2 (6 l /2) and A 10 = -sin 2 (6»/2). The general solution follows 
as 

Poo = 

— £00,11 [£01,01 ~~ £01,10 + ^ 2 (A)] + 2£ o,oi£oi, 11 (£01,01 — £01,10) 



(£00,00 — £00,11) [£qi,oi ~ £01,10 + ^ 2 (^)] — 2£ o,oi (£01,00 — £01,11) (£01,01 — £ 



01,10; 



Rep, 



(£01,01 — £01,10) [£01,11 + (£01,00 ~ £oi,n)P' 



00J 



01 72 



£01,01 £01,10 + ^ 2 (A) 

Impoi = r Repoi , (42) 

^01,01 — ^01,10 

where fi(A) = e — e 1 . 

In the following, we calculate the amplitude A according to Eq. (j2HJ) to zero-th order 
in e. In Fig. El we show the nonlinear response for the parameter set (in dimensionless 
units) / = 10 -5 and v = 10~ 3 . Moreover, the one-phonon resonance condition reads 
u ex = U\ + v. The transition from the resonant to antiresonant behavior depends on 
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the ratio 7/^10 = 7/(2/). For the case of stronger damping 7/(2/) = 10, we find 
that the response shows a resonant behavior with a Lorentzian form similar to the 
response of a damped linear oscillator. In fact, the corresponding standard classical 
result is also shown in Fig. El (black dashed line). The only effect of the nonlinearity to 
lowest order perturbation theory is to shift the resonance frequency by the nonlinearity 
parameter v. The resonant behavior turns into an antiresonant one if the damping 
constant is decreased to smaller values. A cusp-like line profile arises in the limit of very 
weak damping when the damping strength is smaller than the minimal splitting, i.e., 
7/(2/) <C 1. Then, the response follows from the master equation ()40|) as 



.4 



x 



V2 



e e 

sin - cos - 
2 2 



sin 



cos 



sin 4 I + cos 4 I 



(43) 



2 1 2 

This antiresonance lineshape is also shown in Fig. |2] (see dotted-dashed line). At 
resonance A = 0, we have an equal population of the quasienergy states: p o — P11 = 1/2 
and both add up to a vanishing oscillation amplitude A since A Q0 = —An. Note that we 
show also the solution from the exact master equation containing all orders in e, for the 
case 7/(2/) = 0.5 and s = 1 (blue dashed line in Fig. |2J), in order to verify the validity 
of our perturbative treatment. 



6.2. Multiphonon resonance vs. antiresonance 

In this subsection we want to investigate the multiple multiphonon resonances N > 1. 
In order to illustrate the physics, we start with the simplest case at resonance and within 
the secular approximation. 



6.2.1. Secular approximation at resonance The zero-th order quasienergy solutions are 
given in terms of the eigenstates of the harmonic oscillator in Eq. (jlfij) with 9 = 7r/4. 
Then, \n) and \N — n) (n < N/2) form a pair of quasienergy solutions. For N odd, there 
are (N + l)/2 pairs. For N even, there are N/2 pairs whereas the state \4>n/2) — 1^/2) 
remains sole. Within the secular approximation, we can describe the dynamics in terms 
of the kinetic equation f40|) . Plugging Eq. (|T6|) into the expression for the transition 
rates in Eq. (|4~T]) . we find that most of the transition rates between two different states 
belonging to two different pairs are zero, except for 

^nn, ra+ln+1 £"nn,N— n— liV— n— 1 ^~ , N—nN—n,n+ln+l 

7s 

£n+ln+l,nn £"n+ln+l,N— nN— n £iV— n— IN— n— l,nn 

7s 

= ^■N-n-lN-n-l,N-nN-n = — n ) , 

£n/2N/2,n/2±in/2±i = -j(N + 2) (for N even) , 

Cn/2±in/2±i,n/2N/2 = -jN (for Af even). (44) 
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The transition rates between states belonging to the same pair are zero with the 

exception £(jv-i)/2(jv-i)/2,(jv+i)/2(jv+i)/2 = ls{N + l)/8. 

The dynamics can be illustrated with a simple analogy to a double-well potential. 
Each partner of the pair \<f> n ) and |0jv_„) of the quasienergy states consists of a 
superposition of two harmonic oscillator states |n) and \ N—n) which are the approximate 
eigenstates of the static anharmonic potential in the regime of weak nonlinearity. In 
our simple picture, \n) and \N — n) should be identified with two localized states in 
the two wells of the quasienergy potential, see Fig. El for illustration. Note that a 
quasipotential can be obtained by writing the RWA Hamiltonian in terms of the two 
canonically conjugated variables X and V [H]. The right/left well should be identified 
with the internal/external part of the quasieneryg surface shown in Ref. [4*T| . 

In the figure, we have chosen N = 8. Within our analogy, the states 
|0), |JV/2-l) are localized in one (here, the left) well, while \N), \N-1),..., \N/2 + 

1) are localized in the other well (here, the right). The fact that the true quasienergy 
states are superpositions of the two localized states is illustrated by a horizontal arrow 
representing tunneling. 

From Eq. (|44j) follows that a bath-induced transition is only possible between states 
belonging to two different neighboring pairs. As discussed after Eq. ([41 p. the only 
contribution to the transition rates come from nearby HO states. In our case, we 
consider only spontaneous emission which corresponds to intrawell transitions induced 
by the bath. This is shown schematically in Fig. |3] by the vertical arrows with their 
thickness being proportional to the transition rates. We emphasize that the bath- 
induced transitions occur towards lower lying HO states. Consequently, in our picture, 
spontaneous decay happens in the left well downwards but in the right well upwards. 

The driving field excites the transition from |0) to \N) while the bath generates 
transitions between HO states towards lower energies according to \N) — > \N — 1) — > 
... — > |0) when only spontaneous emission is considered. 

As a consequence, the ratio of the occupation numbers of two states belonging to 
two neighboring pairs is simply given by the ratio of the corresponding transition rates 
according to 



Hence, the unpaired state \(j>N/2) (for N even) or the states \4>(n-i)/2) and |</>(jv+i)/2) (for 
N odd) are the states with the largest occupation probability. By iteration, one finds 



Qnn QN—nN—n i 

Qnn -^nn/re+ln+l 



71+ 1 



(45) 



Qn+ln+1 '~"n,-{-Xn-\-X,nn 



N -n 




i 



0.5,0.37,0.31,0.27,... 



for iV = 2,4,6,8,... , (46) 
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and 



(N-l)/2 n _i 



Q { N T l)/2= 2 + 2 n 



n=l k=0 



N — I — 2k 

N + 3 + 2k 



■■ 0.37, 0.31, 0.27, 0.25,... 
for N = 3, 5, 7, 9, . . . (47) 



Density matrix around the resonance So far, we have discussed the dynamics 
exactly at resonance. Next, we consider the situation around the resonance and for an 
increased coupling to the bath. Therefore, we compute the stationary solution using 
the conditions in Eq. and the general leading order solution for the quasienergy 
states given in Eq. f|16j) . The expressions for the rates which are modified compared to 
before follow straightforwardly and are given in the Appendix. Similarly, the only three 
equations which change compared to the previous situation are also presented there. 
These equations can be straightforwardly solved by 

2 



£?oo 



Qnn 



Qno 



1 9 9 N 

— cot - H 

N 2 2 



i 2 e n 

iv tan 2 + y 



7s 



fi(A) 

7s 
fi(A) 



2$ 
COS - 



' • 20 
Sill - 

2 



. 
sin - cos - 
2 2 



iV 



7s 



7s 



1 9 9 1 9 9 
1 + - tan 2 - + - cot 2 - 

2 2 2 2 



1 2 6» 1 2 6> 
1 + - tan 2 - + - cot 2 - 

2 2 2 2 

1 2 # 1 2 # 

1 + - tan 2 - + - cot 2 - 

2 2 2 2 



£11 



£11 



(4f 



(49) 



|0)(0|, (50) 



2 fi(A) y n(A) 

Away from the resonance (\9\ <C 1), the density matrix follows as 
^|0 o )(0o|~|O)(O|. 

In the limit of strong coupling (j s ^> fi(A)), one finds 

9 9 9 9 

q ~ cos 2 -|0o)(0o| + sin 2 -\(Pn)((Pn\ + sin - cos -(|0o)(<M + |0jv)(0o|) 

for any ^. This nicely illustrates that when the coupling to the bath is strong enough, 
the possibility of resonant tunneling between |0) and \N) is destroyed and a trivial 
asymptotic state results. This is true even if tunneling transitions between the other 
states are possible. Moreover, this also shows that moving away from resonance also 
suppresses multiphonon tunneling transitions. In other words, the only requirement 
for the multiphonon transition to occur in the stationary limit is the possibility of the 
tunneling transition |0) — > \N). 



6.2.3. Lineshape around the resonance Within our partial secular approximation, the 
lineshape of the oscillator nonlinear response given in Eq. (|3fjjl reduces to 

A = \/2x \ g aa A aa + g 0N A 0N + g N0 A N0 \ . (51) 

a/3 

The leading order is given by the zero-th order expression for g and the first-order 
expressions for A aa ,A N0 and A 0N . In order to compute these matrix elements, we 
determine the first order eigenvectors using Van Vleck perturbation theory according to 

|0o)i = e ie5l |0 o )o, (52) 
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where S\ is the first order component in the expansion of S with respect to e given in 
Eq. (jllj) . The matrix elements of its off-diagonal blocks are given by 

(pi\Si\f3) = i— — . (53) 

Here, E a are the eigenenergies of the unperturbed Hamiltonian Hq given in Eq. (fTUj) . 
This yields for N = 2 

(6 6\ ( 6 6 

1 - 2^2sin - cos - j , A 22 = 3e ( 1 + 2^ sin - cos - 

6 6 
A 02 = 6V2e cos 2 - A 20 = -6V2e sin 2 - , A u = -9e . (54) 

The corresponding result for the nonlinear response for iV = 2 is shown in Fig. 0] for 
the case v = and / = 10~ 4 for different values of 7 s /fi2o- F° r strong damping 
7s/^2o — 5, the resonance is washed out almost completely. Decreasing damping, 
a resonant lineshape appears whose maximum is shifted compared to the resonance 
condition u cx = ui + 3b , /2. Note that the dashed line refers to the result which includes 
all orders in e and which follows from the numerical solution of the master equation for 
an Ohmic bath at temperature T = 0.1T . The picture which arises for the behavior 
is the following: For weak damping (7^ <C ^20); the equilibrium state is a statistical 
mixture of quasienergy states. At resonance, the most populated state is \4>i) which 
oscillates with a phase difference of — ir in comparison with the driving. This is due to 
the negative sign of An in Eq. (|54jl . Hence, at resonance the overall oscillation of the 
observable occurs with a phase difference of ip = —it. Far away from resonance, the 
most populated state is \<po), see Eq. (f4*H|) . which oscillates in phase with the driving. 
Thus, the overall oscillation occurs in phase, i.e., <p = 0. If no off-diagonal element of the 
density matrix is populated (which is the case for weak damping), the overall phase is 
either ip = or ip = —71. Hence, increasing the distance from resonance, the amplitude 
A has to go through zero yielding a cusp-like line-shape. This implies the existence of 
a maximum in the response. For slightly larger damping, the finite population of the 
off-diagonal elements leads to a smearing of the cusp. For larger damping, the resonance 
is washed out completely, as has been already discussed, see Eq. (JHDJ). In this regime, 
the oscillation is in phase with the driving. By decreasing the damping, the population 
of the out-of-phase state starts to increase near the resonance resulting in a reduction 
of the in-phase-phase oscillation and thus producing a minimum of the response. 

This mechanism is effective for a broad range of parameters including larger N, 
larger e, and larger temperature T, as also shown numerically in Refs. [SU EHj- Note 
that a calculation up to first order in e for the density matrix is required for iV odd, since 
the matrix elements j4(jv±i)/2(jv±i)/2 have a zero-th order term, in order that the overall 
result for A is again of first order in e. Since one obtains more complicated expressions 
than before, we omit to present them in their full lengths. In Fig. we show the 
behavior for N = 3 for various damping constants 7^/^30 for the case / = 0.5 x 10~ 4 
and v = 1CT 3 . For a large value for 7 s /^3o, the resonance is washed out completely. 
When the damping is decreased, a dip appears which corresponds to an antiresonance. 



Nonlinear response of a driven vibrating nanobeam in the quantum regime 19 

Decreasing the damping further, the antiresonance turns into a clear resonance. This 
behavior is opposite to the case N = 1 as discussed above, but similar to the case N = 2. 

7. Conclusions 

We have studied the nonlinear response of a vibrating nanomechanical beam to a time- 
dependent periodic driving. Thereby, a static longitudinal compression force is included 
and the system is investigated close to the Euler buckling instability. There, the 
fundamental transverse mode dominates the dynamics whose amplitude can be described 
by an effective single-particle Hamiltonian with a periodically driven anharmonic 
potential with a quartic nonlinearity. Damping is modeled phenomenologically by a 
bath of harmonic oscillators. We allow for an Ohmic as well as for a super-Ohmic 
spectral density and have considered the regime of weak system-bath coupling. In this 
regime, the dynamics is captured by a Born-Markovian master equation formulated 
in the frame which rotates with the driving frequency. The pure driven Hamiltonian 
shows avoided level crossings of the quasienergies which correspond to multiphonon 
transitions in the resonator. In fact, a transition between a resonant and an antiresonant 
behavior at the avoided level crossings has been found which depends on the coupling to 
the bath. Concentrating to driving frequencies around the avoided level crossings, the 
dynamics can be simplified considerably by restricting to a few quasienergy levels. In 
order to illustrate the basic principles governing the resonance-antiresonance transition, 
we investigate the perturbative regime of weak nonlinearity and weak driving strength. 
Then, Van Vleck perturbation theory allows to calculate the quasienergies and the 
quasienergy states and an analytic solution of the master equation becomes possible 
yielding directly the nonlinear response. For the one-phonon case, we find a simple 
expression for the nonlinear response which displays a Lorentzian resonant behavior 
for strong damping. Reducing the damping strength, an antiresonance arises. For the 
multiphonon transitions, first an antiresonance arises when damping is reduced. For 
even smaller values of the damping constant, the antiresonance turns into a resonant 
peak. This is due to a subtle interplay of varying populations of quasienergy states 
which is affected by the bath. 

Finally, a comment on the observability of this effect is in order. The amplitude A 
measuring the nonlinear response is of the order of the oscillator length scale xq. This 
makes it challenging to measure the effect directly since the deterministic vibrations are 
on the same length scale like the quantum fluctuations. In turn, more subtle detection 
strategies have to be worked out, for instance, the capacitive coupling of the resonator 
to a Cooper pair box jHUIEE] or single electron transistors jMllSZl, the use of squeezed 
states in this setup [22], the use of a second nanoresonator as a transducer for the phonon 
number in the first one [37] , or the coherent signal amplification by stochastic resonance 
[HE]. In any case, the experimental confirmation of the theoretically predicted effects 
remains to be provided. 
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Appendix: Density matrix around the multiphonon resonance 

For completeness, we present in this Appendix the calculation of the density matrix 
around the multiphonon resonance which is required for Section 16.2.21 The expressions 
for the rates which are modified compared to before are readily calculated to be 

7s 2 

oo,jv-uv-i — -n cos - , 

7s • 2 6 
NN,N-1N-1 — ~2 Sm 2 ' 



£00,11 
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£aw,h 
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£11,00 
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£ll,AW 
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£oo,oiv 
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(55) 



Similarly, there are only three equations which change compared to the previous 
situation. They read 

9 6 6 6 

= - iVsin 2 -£oo + cos 2 -g u + A cos - sin -g m , 







6 6 6 6 

N cos 2 — Qnn + sin 2 ~g n + N cos - sin -g N0 



= -iQ(A)g N0 + 



Is 



6 6 

-Ng N0 + cos - sin -(iV^oo + Ng NN + 2g n ) 



(56) 



with the quasienergy level splitting 

n(A) =e N -e = -sgn(A) 



+^ 2 



no 



1/2 
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Figure 1. Typical quasienergy spectrum e a for increasing driving frequency w ex for 
the case v — 10~ 3 and / = 10~ 4 . The vertical dashed lines indicate the multiple 
avoided level crossings for a fixed driving frequency. 
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Figure 2. Nonlinear response of the nanoresonator at the one-phonon resonance 
N = 1, where cu cx = uii+v for the parameters / = 10 -5 and v = 10~ 3 (in dimensionless 
units). The transition from a resonant behavior for large damping (7/(2/) = 10) to an 
antiresonant behavior at small damping (7/(2/) -C l)is clearly visible. The resonant 
line shape is a Lorentzian and coincides with the linear response of a harmonic oscillator 
at frequency u)\ + v (see black dashed line for 7/(2/) = 10). Also shown is the limit 
of 7/(2/) -c 1 (black dotted-dashed line) yielding a cusp-like lineshape. Note that we 
also depict the solution from the exact master equation for the case 7/(2/) = 0.5 (blue 
dashed line). 
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Figure 3. Schematic view of the quasipotential and localized states in the rotating 
frame for the case N — 8. Shown are the pairs of harmonic oscillator states consisting 
of \n) and | N — n) each of which is localized in one of the two wells. The corresponding 
quasienergy states \cf> n } and \(j>N- n ) are a superposition of the two localized states, 
see text. The horizontal arrows indicate the multiphonon transitions between the two 
quasienergy states. The vertical arrows mark the bath-induced transitions with their 
thickness being proportional to the transition rate. 
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Figure 4. Nonlinear response at the two-phonon resonance N = 2, where u> cx = 
lui + 3i//2 for the parameters / = 10~ 4 and v = 10~ 3 (in dimensionless units) for 
different values of the damping. 
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Figure 5. Nonlinear response at the three-phonon resonance N — 3, where cu cx — 
u>i + 2v for the parameters / = 0.5 x 10 -4 and v — 10~ 3 (in dimensionless units) for 
different values of the damping. 



